Appendix 3: Additional models results
#Loading data
load("scripts/datasets.Rdata")
load("models_outputs/models_combined.Rdata") # saved models results combined traits
load("models_outputs/models_traits.Rdata") # saved models results separate traits#Loading auxiliary function for R2 calculations
source("scripts/r2_auxiliary_function.R")# variance inflation factor function
vif.mer <- function (fit) {
## adapted from rms::vif
v <- vcov(fit)
nam <- names(fixef(fit))
## exclude intercepts
ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
if (ns > 0) {
v <- v[-(1:ns), -(1:ns), drop = FALSE]
nam <- nam[-(1:ns)]
}
d <- diag(v)^0.5
v <- diag(solve(v/(d %o% d)))
names(v) <- nam
v
}1. Models with traits
Specification of the models: We used lme4 package to perform a GLMM with binomial (proportion) distribution. An example of the code for each dataset are as follows:
mhigh.spe <- glmer(cbind(occor, n.visit-occor) ~
forest_site400*lbody_size +
forest_site400*nest +
forest_site400*diet +
forest_site400*lower_stratum +
forest_land*lbody_size +
forest_land*nest +
forest_land*diet +
forest_land*lower_stratum +
(forest_site400 + forest_land|sp) +
(1|landscape:sp) + (1|site:sp) +
(lbody_size + nest + diet + lower_stratum|landscape) +
(lbody_size + nest + diet + lower_stratum|site),
family=binomial, data=high.spe,
nAGQ = 1, control = glmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 500000)))
We ran separate models for each assemblage and trait. Afterwards, we ran one model with the combination of the traits body mass, diet, nest type and % of lower strata use. Table S3. shows the marginal R2 of all models terms.
traits = c("lbody_size", "nest", "diet", "lower_stratum")
rhigh.spe <- rbind(r2.calc(mhigh.spe, atrib = traits) ,
r2.calc(bodyhigh.spe, atrib="lbody_size"),
r2.calc(nesthigh.spe, atrib="nest"),
r2.calc(diethigh.spe, atrib="diet"),
r2.calc(frugihigh.spe, atrib="frugivory"),
r2.calc(insehigh.spe, atrib="insectivory"),
r2.calc(low_strathigh.spe, atrib = "lower_stratum"),
r2.calc(stratumhigh.spe, atrib="stratum"))
rhigh.spe <- cbind(model = rep(c("Combined", "body mass", "nest type", "main diet",
"% frugivory", "% insetivory","% lower strata",
"foraging stratum")),
rhigh.spe) %>%
mutate_if(is.numeric, function(x) 100*x)
colnames(rhigh.spe) <- c("Model", "Total", "trait*env", "env|sp", "lands:sp",
"site:sp", "trait|lands", "trait|site")
rlow.spe <- rbind(r2.calc(mlow.spe, atrib = traits) ,
r2.calc(bodylow.spe, atrib="lbody_size"),
r2.calc(nestlow.spe, atrib="nest"),
r2.calc(dietlow.spe, atrib="diet"),
r2.calc(frugilow.spe, atrib="frugivory"),
r2.calc(inselow.spe, atrib="insectivory"),
r2.calc(low_stratlow.spe, atrib = "lower_stratum"),
r2.calc(stratumlow.spe, atrib="stratum"))
rlow.spe <- cbind(model = rep(c("Combined", "body mass", "nest type", "main diet",
"% frugivory", "% insetivory","% lower strata",
"foraging stratum")),
rlow.spe) %>%
mutate_if(is.numeric, function(x) 100*x)
colnames(rlow.spe) <- c("Model", "Total", "trait*env", "env|sp", "lands:sp",
"site:sp", "trait|lands", "trait|site")
rhigh.gen <- rbind(r2.calc(mhigh.gen, atrib = traits) ,
r2.calc(bodyhigh.gen, atrib="lbody_size"),
r2.calc(nesthigh.gen, atrib="nest"),
r2.calc(diethigh.gen, atrib="diet"),
r2.calc(frugihigh.gen, atrib="frugivory"),
r2.calc(insehigh.gen, atrib="insectivory"),
r2.calc(low_strathigh.gen, atrib = "lower_stratum"),
r2.calc(stratumhigh.gen, atrib="stratum"))
rhigh.gen <- cbind(model = rep(c("Combined", "body mass", "nest type", "main diet",
"% frugivory", "% insetivory","% lower strata",
"foraging stratum")),
rhigh.gen) %>%
mutate_if(is.numeric, function(x) 100*x)
colnames(rhigh.gen) <- c("Model", "Total", "trait*env", "env|sp", "lands:sp",
"site:sp", "trait|lands", "trait|site")
rlow.gen <- rbind(r2.calc(mlow.gen, atrib = traits) ,
r2.calc(bodylow.gen, atrib="lbody_size"),
r2.calc(nestlow.gen, atrib="nest"),
r2.calc(dietlow.gen, atrib="diet"),
r2.calc(frugilow.gen, atrib="frugivory"),
r2.calc(inselow.gen, atrib="insectivory"),
r2.calc(low_stratlow.gen, atrib = "lower_stratum"),
r2.calc(stratumlow.gen, atrib="stratum"))
rlow.gen <- cbind(model = rep(c("Combined", "body mass", "nest type", "main diet",
"% frugivory", "% insetivory","% lower strata",
"foraging stratum")),
rlow.gen) %>%
mutate_if(is.numeric, function(x) 100*x)
colnames(rlow.gen) <- c("Model", "Total", "trait*env", "env|sp", "lands:sp",
"site:sp", "trait|lands", "trait|site")todos <- bind_rows(list(rhigh.spe, rlow.spe, rhigh.gen, rlow.gen))
kable(todos, "latex", booktabs=T, caption= "Overall and marginal R-squared of trait models in each dataset. For the marginal R-squared terms see Table 2 (main text).") %>%
group_rows("Specialists", 1, 16) %>%
group_rows(" High quality", 1,8) %>%
group_rows(" Low quality", 9,16) %>%
group_rows("Generalists", 17, 32) %>%
group_rows(" High quality", 17,24) %>%
group_rows(" Low quality", 25,32)2. Models coeficients
Tables S3., S3., S3., and S3. show the coefficients for each model.
tidy(mhigh.spe, effects = "fixed") %>%
kable(digits=2,"latex", booktabs=T, caption = "Fixed effects oefficients for the model of specialists in high-quality matrix landscapes.") tidy(mlow.spe, effects = "fixed") %>%
kable(digits=2,"latex", booktabs=T, caption = "Fixed effects oefficients for the model of specialists in low-quality matrix landscapes.") tidy(mhigh.gen, effects = "fixed") %>%
kable(digits=2, "latex", booktabs=T, caption = "Fixed effects oefficients for the model of generalists in high-quality matrix landscapes.") tidy(mlow.gen, effects = "fixed") %>%
kable(digits=2,"latex", booktabs=T, caption = "Fixed effects oefficients for the model of generalists in low-quality matrix landscapes.")3. Models diagnostic
Variance Inflation Factor of the model parameters for each dataset in Table S3..
vifhigh.spe <- data.frame(termo = names(vif.mer(mhigh.spe2)),
high.spe =vif.mer(mhigh.spe2))
viflow.spe <- data.frame(termo = names(vif.mer(mlow.spe2)),
low.spe =vif.mer(mlow.spe2))
viflow.gen <- data.frame(termo = names(vif.mer(mhigh.gen2)),
low.gen =vif.mer(mhigh.gen2))
viflow.gen <- data.frame(termo = names(vif.mer(mlow.gen2)),
low.gen =vif.mer(mlow.gen2))
varif <- vifhigh.spe %>% full_join(viflow.spe, "termo") %>%
full_join(viflow.gen, "termo") %>% full_join(viflow.gen, "termo")
varif$termo <- c("forest.local", "body_mass", "nest_closed", "nest_open_semi",
"diet_insectivorous", "diet_onivorous", "lower_strata",
"diet_granivorous", "forest.landscape", "diet_nectarivorous")
colnames(varif) <-c("parameter","Coffee", "lowture", "Coffee", "lowture")
kable(varif, "latex", booktabs=T, caption= "Variance Inflation Factor index for combined traits models in each dataset.", digits=2) %>%
add_header_above(c(" " = 1, "Specialists" = 2, "Generalists"= 2))Example of the residual diagnostic of the model with the combined traits (main diet, body mass, nest type and % of lower strata use) for the forest specialists in high-quality matrix landscapes. The models’ diagnostics for the other assemblages were all similar and can be checked in this Rmd file.
Residual correlations among species and sites
Below we present the Kendall correlations for the residuals among species and sites for the models using the predictions for site:sp random effect (Observation Level Random Effect). For the residual correlations we followed the code provided by Miller, Damschen & Ives (2018).
nsites = 40
nspp = length(unique(high.spe$sp))
dat <- high.spe
dat$resid <- as.matrix(ranef(mhigh.spe)$`site:sp`)
X <- matrix(dat$resid, nrow=nsites, ncol=nspp, byrow=F)
rownames(X) <- unique(dat$ponto)
colnames(X) <- unique(dat$sp)
# species correlations
corrs.sp <- cor(X, method="kendall")
for(i in 1:nspp) corrs.sp[i,i] <- NA
# site correlations
corrs.site <- cor(t(X), method="kendall")
for(i in 1:nsites) corrs.site[i,i] <- NARange of species correlations: -0.4, 0.43. Range of sites correlations: -0.3, 0.27.
corrplot(corrs.sp, method="square", type="upper",diag=F,tl.cex=0.5,cl.cex=0.5)Species residual Kendall correlations for the specialist species in high-quality matrix landscapes.
corrplot(corrs.site, method="square", type="upper",diag=F,tl.cex=0.5,cl.cex=0.5)Sites residual Kendall correlations for the specialist species in high-quality matrixlandscapes.
x.sp <- matrix(corrs.sp, ncol=1)
x.site <- matrix(corrs.site, ncol=1)
x.sp <- x.sp[!is.na(x.sp)]
x.site <- x.site[!is.na(x.site)]
# Figure histogram
par(mfrow=c(1,2), mai=c(1,1,.5,.2))
hist(corrs.sp, main="", probability=T, breaks=.1*(-10:10), xlab="Correlation coefficients", ylab="Probability", ylim=c(0,4), cex.lab=1.2)
mtext("Species", side=3,cex=1.2)
hist(corrs.site, main="", probability=T, breaks=.1*(-10:10), xlab="Correlation coefficients", ylab="Probability", ylim=c(0,4), cex.lab=1.2)
mtext("Sites", side=3, cex=1.2)Histograms of the residual Kendall correlations for the specialists species in high-quality matrix landscapes.
Residual diagnostic
We used DHARMa package (Hartig (2018)) for the diagnostic of quantile residuals.
resce <- simulateResiduals(mhigh.spe, n=1000, ref.form=~0)
plot(resce)Residuals against predictors:
par(mfrow=c(4,2))
plotResiduals(resce, high.spe$forest_site400, xlab="Local forest cover 400m")
plotResiduals(resce, high.spe$forest_land, xlab="Landscape forest")
plotResiduals(resce, high.spe$lbody_size, xlab="Body size")
plotResiduals(resce, high.spe$nest, xlab="Nest type")
plotResiduals(resce, high.spe$diet, xlab="Main diet")
plotResiduals(resce, high.spe$low_strat, xlab="% lower strata use")Predictions for each species local forest cover
Landscape forest cover was fixed in 30%.
#landscape cover 30%
fland = (30 - mean(high.spe$forest_landOrig))/
sd(high.spe$forest_landOrig)
# local forest cover
f400st.high <- seq(min(high.spe$forest_site400), max(high.spe$forest_site400),
length.out=5)
f400Orig.high <- f400st.high * sd(high.spe$forest_site400Orig) +
mean(high.spe$forest_site400Orig)
local.high <- data.frame(forest_site400=f400st.high,
forest_site400Orig = f400Orig.high)
f400st.low <- seq(min(low.spe$forest_site400), max(low.spe$forest_site400),
length.out=5)
f400Orig.low <-f400st.low * sd(low.spe$forest_site400Orig) +
mean(low.spe$forest_site400Orig)
local.low <- data.frame(forest_site400=f400st.low,
forest_site400Orig = f400Orig.low)# high quality specialists
sp.tra <- high.spe %>% select(sp, lbody_size, nest, diet, lower_stratum) %>%
distinct()
shigh.spe <- expand.grid(sp = unique(high.spe$sp),
forest_land = fland,
forest_site400 =f400st.high,
landscape = "P02",
site = "P02.P00") %>%
left_join(sp.tra, "sp") %>%
left_join(local.high, "forest_site400")
shigh.spe$pretes <- predict(mhigh.spe, newdata=shigh.spe, type="response",
re.form=NULL)
# low quality specialists
sp.tra2 <- low.spe %>% select(sp, lbody_size, nest, diet, lower_stratum) %>%
distinct()
slow.spe <- expand.grid(sp = unique(low.spe$sp),
forest_land = fland,
forest_site400 =f400st.low,
#forest_site400Orig = f400Orig.low,
landscape = "148",
site = "148.P10") %>%
left_join(sp.tra2, "sp")%>%
left_join(local.low, "forest_site400")
slow.spe$pretes <- predict(mlow.spe, newdata=slow.spe, type="response",
re.form=NULL)bind_rows(list(High=shigh.spe,Low=slow.spe), .id="Matrix quality") %>%
mutate(forest_site400Orig = round(forest_site400Orig)) %>%
ggplot(aes(x=forest_site400Orig, y=pretes, col=`Matrix quality` )) +
geom_line() +
facet_wrap(~sp, ncol=8) +
theme_cowplot() +
theme(legend.position = "top",
strip.text.x = element_text(hjust=0, size=10))+
scale_y_continuous(name= "Occurrence probability") +
scale_x_continuous(name="Local forest cover (%)") +
ggtitle("Forest specialists")# high quality generalists
sp.tra <- high.gen %>% select(sp, lbody_size, nest, diet, lower_stratum) %>%
distinct()
shigh.gen <- expand.grid(sp = unique(high.gen$sp),
forest_land = fland,
forest_site400 =f400st.high,
landscape = "P02",
site = "P02.P00") %>%
left_join(sp.tra, "sp") %>%
left_join(local.high, "forest_site400")
shigh.gen$pretes <- predict(mhigh.gen, newdata=shigh.gen, type="response",
re.form=NULL)
# low quality generalists
sp.tra2 <- low.gen %>% select(sp, lbody_size, nest, diet, lower_stratum) %>%
distinct()
slow.gen <- expand.grid(sp = unique(low.gen$sp),
forest_land = fland,
forest_site400 =f400st.low,
#forest_site400Orig = f400Orig.low,
landscape = "148",
site = "148.P10") %>%
left_join(sp.tra2, "sp")%>%
left_join(local.low, "forest_site400")
slow.gen$pretes <- predict(mlow.gen, newdata=slow.gen, type="response",
re.form=NULL)bind_rows(list(High=shigh.gen,Low=slow.gen), .id="Matrix quality") %>%
mutate(forest_site400Orig = round(forest_site400Orig)) %>%
ggplot(aes(x=forest_site400Orig, y=pretes, col=`Matrix quality` )) +
geom_line() +
facet_wrap(~sp, ncol=8) +
theme_cowplot() +
theme(legend.position = "top",
strip.text.x = element_text(hjust=0, size=10))+
scale_y_continuous(name= "Occurrence probability") +
scale_x_continuous(name="Local forest cover (%)") +
ggtitle("Forest generalists")